!> \brief
!! This is a testcase for the Stokes problem based on the squeezing of
!! a fluid between two plates. This test case has an analytic
!! solution.
!!
!! \n
!!
!! Th computational domain is
!! \f$\Omega=[-A,A]\times[-B,B]\times[-s,s]\f$, with the
!! two-dimensional case obtained for \f$B\to\infty\f$. The viscosity
!! is \f$\mu=1\f$ and there is no forcing term. The analytic solution
!! is
!! \f{displaymath}{
!!  \underline{u} = \left[\begin{array}{c}
!!   \alpha x \left(1-\frac{z^2}{s^2}\right) \\[3mm]
!!   \beta y \left(1-\frac{z^2}{s^2}\right) \\[3mm]
!!   -\frac{W_0}{2s}z\left(3-\frac{z^2}{s^2}\right)
!!  \end{array}\right], \qquad
!!  p = \frac{1}{2}\frac{W_0}{s^3}\left(
!!   \frac{B^2}{A^2+B^2}(A^2-3x^2)
!!  +\frac{A^2}{A^2+B^2}(B^2-3y^2)
!!  -(s^2-3z^2)\right),
!! \f}
!! where \f$W_0\in\mathbb{R}\f$ and
!! \f{displaymath}{
!!  \alpha = \frac{3}{2}\frac{B^2}{A^2+B^2}\frac{W_0}{s}, \qquad
!!  \beta = \frac{3}{2}\frac{A^2}{A^2+B^2}\frac{W_0}{s}.
!! \f}
!!
!! On the \f$z=\pm s\f$ boundaries a Dirichlet boundary condition
!! should be enforced, while on the remaining boundaries either
!! Dirichlet or Neumann boundary conditions can be specified; for
!! specifying Neumann boundary conditions, the stress tensor is
!! \f{displaymath}{
!!  \underline{T} =
!!  \mu\left(\nabla\underline{u}+(\nabla\underline{u})^T\right)
!!  -p\underline{I} =
!!   \left[\begin{array}{ccc}
!!    2\alpha\xi & 0 & -\frac{2\alpha}{s^2}xz \\[3mm]
!!    0 & 2\beta\xi  & -\frac{2\beta}{s^2}yz \\[3mm]
!!    -\frac{2\alpha}{s^2}xz & -\frac{2\beta}{s^2}yz  &
!!    -3\frac{W_0}{s}\xi
!!  \end{array}\right] - p\underline{I},
!! \f}
!! with \f$\xi=1-\frac{z^2}{s^2}\f$. Finally, we notice that the curl
!! of the velocity is
!! \f{displaymath}{
!!  \nabla\times\underline{u} = \left[\begin{array}{c}
!!   \frac{2\beta}{s^2}yz \\[3mm]
!!   -\frac{2\alpha}{s^2}xz \\[3mm]
!!    0
!!  \end{array}\right].
!! \f}
!<----------------------------------------------------------------------
module mod_squeeze_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_squeeze_test_constructor, &
   mod_squeeze_test_destructor,  &
   mod_squeeze_test_initialized, &
   test_name, &
   test_description,&
   coeff_visc,&
   coeff_f,   &
   coeff_dir, &
   coeff_u,   &
   coeff_gradu,&
   coeff_p,   &
   coeff_gradp

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter ::    &
   test_name = "squeeze"
 ! private members
 real(wp), parameter :: &
   a = 2.0_wp, &
   b = 3.0_wp, &
   s = 1.0_wp, &
   w0 = 1.0_wp

! Module variables

 ! public members
 character(len=100), protected ::    &
   test_description(2)
 logical, protected ::               &
   mod_squeeze_test_initialized = .false.
 ! private members
 integer :: d ! dimension
 real(wp) :: alpha, beta
 character(len=*), parameter :: &
   this_mod_name = 'mod_squeeze_test'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_squeeze_test_constructor(id)
  integer, intent(in) :: id

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_squeeze_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   d = id

   write(test_description(1),'(a,i1)') &
     'Squeeze test case, dimension d = ',d
   testdim: select case(d)
    case(2)
     alpha = 3.0_wp/2.0_wp * w0/s
     write(test_description(2),'(a,e9.3,a,a,a,e9.3,a,e9.3,a)') &
       '  parameters: A = ',a,', B = ','Inf',', s = ',s,', W0 = ',w0,'.'
    case(3)
     alpha = 3.0_wp/2.0_wp * b**2/(a**2+b**2) * w0/s
     beta  = 3.0_wp/2.0_wp * a**2/(a**2+b**2) * w0/s
     write(test_description(2),'(a,e9.3,a,e9.3,a,e9.3,a,e9.3,a)') &
       '  parameters: A = ',a,', B = ',b,', s = ',s,', W0 = ',w0,'.'
    case default
     call error(this_sub_name,this_mod_name, &
       'This test case is only implemented for dimensions 2 and 3.')
   end select testdim

   mod_squeeze_test_initialized = .true.
 end subroutine mod_squeeze_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_squeeze_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_squeeze_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_squeeze_test_initialized = .false.
 end subroutine mod_squeeze_test_destructor

!-----------------------------------------------------------------------
 
 pure function coeff_visc(x) result(nu_x)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: nu_x(size(x,2))

   nu_x = 1.0_wp
 end function coeff_visc
 
!-----------------------------------------------------------------------
 
 pure function coeff_f(x) result(f)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: f(size(x,1),size(x,2))

   f = 0.0_wp
 end function coeff_f
 
!-----------------------------------------------------------------------

 pure function coeff_dir(x,breg) result(d)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: d(size(x,1),size(x,2))

   d = coeff_u(x)
 end function coeff_dir
 
!-----------------------------------------------------------------------

 pure function coeff_u(x) result(u)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: u(size(x,1),size(x,2))

  real(wp) :: xi(size(x,2))

   if(d.eq.2) then
     u(1,:) = alpha*x(1,:)*(1.0_wp-x(2,:)**2/s**2)
     u(2,:) = -w0/(2.0_wp*s)*x(2,:)*(3.0_wp-x(2,:)**2/s**2)
   else
     xi = 1.0_wp - x(3,:)**2/s**2
     u(1,:) = alpha*x(1,:)*xi
     u(2,:) =  beta*x(2,:)*xi
     u(3,:) = -w0/(2.0_wp*s)*x(3,:)*(3.0_wp-x(3,:)**2/s**2)
   endif

 end function coeff_u
 
!-----------------------------------------------------------------------

 pure function coeff_gradu(x) result(gu)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: gu(size(x,1),size(x,1),size(x,2))

  real(wp) :: xi(size(x,2))

   if(d.eq.2) then
     xi = 1.0_wp - x(2,:)**2/s**2
     ! first row
     gu(1,1,:) = alpha*xi
     gu(1,2,:) = -2.0_wp*alpha*x(1,:)*x(2,:)/s**2
     ! second row
     gu(2,1,:) = 0.0_wp
     gu(2,2,:) = -3.0_wp/2.0_wp*w0/s*xi
   else
     xi = 1.0_wp - x(3,:)**2/s**2
     ! first row
     gu(1,1,:) = alpha*xi
     gu(1,2,:) = 0.0_wp
     gu(1,3,:) = -2.0_wp*alpha*x(1,:)*x(3,:)/s**2
     ! second row
     gu(2,1,:) = 0.0_wp
     gu(2,2,:) = beta*xi
     gu(2,3,:) = -2.0_wp*beta*x(2,:)*x(3,:)/s**2
     ! third row
     gu(3,1,:) = 0.0_wp
     gu(3,2,:) = 0.0_wp
     gu(3,3,:) = -3.0_wp/2.0_wp*w0/s*xi
   endif

 end function coeff_gradu
 
!-----------------------------------------------------------------------

 pure function coeff_p(x) result(p)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: p(size(x,2))

   if(d.eq.2) then
     p = 0.5_wp*w0/s**3*(a**2-3.0_wp*x(1,:)**2-(s**2-3.0_wp*x(2,:)**2))
   else
     p = 0.5_wp*w0/s**3*( (b**2/(a**2+b**2))*(a**2-3.0_wp*x(1,:)**2) &
                         +(a**2/(a**2+b**2))*(b**2-3.0_wp*x(2,:)**2) &
                         -                   (s**2-3.0_wp*x(3,:)**2) )
   endif

 end function coeff_p
 
!-----------------------------------------------------------------------

 pure function coeff_gradp(x) result(gp)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: gp(size(x,1),size(x,2))

   if(d.eq.2) then
     gp(1,:) = -3.0_wp*w0/s**3*x(1,:)
     gp(2,:) =  3.0_wp*w0/s**3*x(2,:)
   else
     gp(1,:) = -3.0_wp*w0/s**3*(b**2/(a**2+b**2))*x(1,:)
     gp(2,:) = -3.0_wp*w0/s**3*(a**2/(a**2+b**2))*x(2,:)
     gp(3,:) =  3.0_wp*w0/s**3*                   x(3,:)
   endif

 end function coeff_gradp
 
!-----------------------------------------------------------------------

end module mod_squeeze_test

